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Abstract—Madagascar’s flora and fauna are remarkable both for their diversity and supraspecific endemism. Moreover, 
many taxa contain large numbers of species with limited distributions. Several hypotheses have been proposed to 
explain this high level of microendemism, including 1) riverine barrier, 2) mountain refuge, and 3) watershed contraction 
hypotheses, the latter 2 of which center on fragmentation due to climatic shifts associated with Pliocene/Pleistocene 
glaciations. The Malagasy leaf chameleon genus Brookesia is a speciose group with a high proportion of microendemic taxa, 
thus making it an excellent candidate to test these vicariance scenarios. We used mitochondrial and nuclear sequence data 
to construct a Brookesia phylogeny, and temporal concordance with Pliocene/Pleistocene speciation scenarios was tested 
by estimating divergence dates using a relaxed-clock Bayesian method. We strongly reject a role for Pliocene/Pleistocene 
climatic fluctuations in species-level diversification of Brookesia. We also used simulations to test the spatial predictions 
of the watershed contraction model in a phylogenetic context, independent of its temporal component, and found no 
statistical support for this model. The riverine barrier model is likewise a qualitatively poor fit to our data, but some rela- 
tionships support a more ancient mountain refuge effect. We assessed support for the 3 hypotheses in a nonphylogenetic 
context by examining altitude and species richness and found a significant positive correlation between these variables. 
This is consistent with a mountain refuge effect but does not support the watershed contraction or riverine barrier models. 
Finally, we find repeated higher level east-west divergence patterns 1) between the 2 sister clades comprising the Brookesia 
minima group and 2) within the clade of larger leaf chameleons, which shows a basal divergence between western and 
eastern/northern sister clades. Our results highlight the central role of phylogeny in any meaningful tests of species-level 
diversification theories. [Biogeography; Chamaeleonidae; phylogeny; Pleistocene glaciation; relaxed clock; speciation; 


Squamata.] 


Once part of the southern supercontinent of Gond- 
wana, Madagascar, is currently situated some 400 km 
off the southeastern shore of Africa and is famous 
for a remarkable flora and fauna that are increasingly 
threatened by loss of habitat and other human-induced 
changes. Madagascar’s last subaerial connections to 
Africa and India date to approximately 160 and 90 mil- 
lion years ago, respectively (Storey et al. 1995), although 
limited dispersal to and from other continents may have 
been possible over land bridges or small sea barriers in 
the Late Cretaceous or Paleocene (Briggs 2003; Noonan 
and Chippindale 2006; Van Bocxlaer et al. 2007; van 
der Meijden et al. 2007). This long history of isolation 
has contributed greatly to the remarkable degree of 
floral and faunal endemism of Madagascar (Goodman 
and Benstead 2003), which amounts to 100% in native 
amphibians and 92% in nonmarine nonavian reptiles 
(Glaw and Vences 2007). 

In terrestrial vertebrates, the majority of species diver- 
sity corresponds to a limited number of endemic clades 
that colonized most of the available bioclimatic regions, 
including the eastern rainforests, central highlands, 
western dry deciduous forests, and southern arid spiny 
forests. There are several sister-taxon pairs that follow 
an east-west pattern and possibly arose by specializa- 
tion to these different climatic conditions (the “eco- 
graphic constraint hypothesis”; see Yoder and Heckman 
2006), from groups such as skinks, snakes, boophine 
treefrogs, and spiders (Nussbaum and Raxworthy 1998; 


Nussbaum et al. 1998; Vences and Glaw 2002, 2003; 
Wood et al. 2007). In other groups, the major evolu- 
tionary splits consistently separate a clade occurring 
roughly in the northern fifth of Madagascar from a more 
southern clade, with examples found in chameleons, 
geckos, skinks, and mouse lemurs (Yoder and Heckman 
2006; Boumans et al. 2007). As a general pattern, the 
spatial distribution of species richness in some higher 
taxonomic groups may have been shaped by a latitudi- 
nal and altitudinal middomain effect (Lees et al. 1999). 
Nonetheless, even within the major bioclimatic regions, 
species-level local endemism is high, and a major goal 
of researchers has been to uncover the causes of this 
pattern (Goodman and Benstead 2003 and references 
therein; Vences et al. 2009). 


Potential Climatic and Geographic Barriers 


In their study of biogeographic patterns in the leaf 
chameleons (genus Brookesia) of northern Madagascar, 
Raxworthy and Nussbaum (1995) recognized 5 rainfor- 
est regions delimited largely by altitudinal differences 
and intervening dry forests and characterized by a high 
degree of endemism in these lizards. Raxworthy and 
Nussbaum (1995) hypothesized that Pleistocene climatic 
fluctuations had caused fragmentation of the rainforest, 
resulting in multiple allopatrically distributed sister- 
taxon pairs (hereafter referred to as the mountain refuge 
[MR] hypothesis). These conclusions were not based 
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on explicit phylogenetic hypotheses; sister-taxon pairs 
were assumed based on general morphological 
similarities. 

Wilmé et al. (2006) analyzed Madagascar’s striking 
microendemicity in the context of watersheds. Using 
museum data for more than 35,000 georeferenced land 
vertebrate specimens, they found that watersheds with 
low-elevation headwaters tended to define centers of 
endemism (COEs), whereas those with connections to 
the 3 highest summits in Madagascar (each at >2000 m) 
tended to contain more widespread species. Wilmé et al. 
(2006) reasoned that during periods of Late Tertiary 
and Pleistocene glacial maxima, aridification caused 
contraction of previously widespread mesic environ- 
ments. Lower elevations were more severely affected 
than higher elevations, leading to fragmentation and 
isolation of low-elevation watersheds. In contrast, areas 
with riverine connections to high-elevation source wa- 
ters were buffered from this effect and thus served as 
“retreat-dispersion watersheds” (RDWs), providing a 
means of refuge and recolonization during dry and wet 
cycles, respectively. The authors presented this model of 
Pliocene—Pleistocene fragmentation and refugia (here- 
after referred to as the watershed contraction [WC] 
hypothesis) as the main generator of Madagascar’s 
famously microendemic biota. The Wilmé et al. (2006) 
study was based entirely on extant species distributions 
coupled with climatological, hydrological, and topo- 
graphical variables and incorporated no phylogenetic 
data. Pearson and Raxworthy (2009) found some sup- 
port for the WC hypothesis in lemurs and day geckos 
of the genus Phelsuma (a climatological gradient effect 
was also cited for these lizards), but once again relied 
on distribution patterns alone. 

Finally, the presence of large rivers flowing from the 
highlands either toward the west or the east has also 
been discussed as a factor influencing species formation 
and microendemism (hereafter referred to as the river- 
ine barrier [RB] hypothesis). In western Madagascar, 
these rivers coincide with the limits of distribution areas 
of species or phylogeographic lineages within species 
of lemurs and are therefore likely to constitute signifi- 
cant barriers to gene flow for these animals (Pastorini 
et al. 2003). A similar situation seems to exist in east- 
ern Madagascar and may especially influence species 
occurring at low elevations where rivers are widest 
(Goodman and Ganzhorn 2004; Louis et al. 2005). For 
a more detailed overview of the processes inherent 
to the MR, WC, and RB hypotheses, see Vences et al. 
(2009). 


Brookesia Leaf Chameleons 


The chameleon genus Brookesia is an excellent group 
to test hypotheses of species-level diversification and 
microendemism on Madagascar. These lizards, which 
appear to form the sister taxon of all other chameleons 
(Rieppel 1987; Townsend and Larson 2002), constitute 
one of the largest Malagasy reptile groups (Raxworthy 
and Nussbaum 1995; Glaw and Vences 2007). Most 


Brookesia species have very small ranges (Raxworthy 
and Nussbaum 1995; Raselimanana and Rakotoma- 
malala 2003), with almost half known essentially from 
single localities (Carpenter and Robson 2005). Brookesia 
can be divided into 3 main groups based on morphol- 
ogy. One is represented by just 2 species (B. nasus and 
B. lolontany) that are highly divergent from other Brooke- 
sia by both molecular (Raxworthy et al. 2002; Townsend 
and Larson 2002) and morphological measures (long 
snouts, laterally compressed bodies; Raxworthy and 
Nussbaum 1995). A second larger group is composed of 
approximately 18 species with more robust bodies and 
blunted snouts. The remaining 6 described species have 
taken the already diminutive body form of Brookesia 
to an extreme, with total lengths of about 28-45 mm, 
making them some of the world’s smallest vertebrate 
species (Glaw and Vences 2007). For ease of reference, 
these 3 groups will hereafter be referred to as the B. nasus 
group, the “typical” Brookesia, and the B. minima group, 
respectively. 


Hypothesis Testing 


In this study, we use Brookesia to test the temporal 
and spatial predictions of 3 species-level diversification 
hypotheses for Madagascar (MR, WC, and RB; Vences 
et al. 2009). The first step in this process is to infer a 
phylogeny of Brookesia to identify statistically supported 
sister-species pairs using DNA sequence data from mul- 
tiple mitochondrial and nuclear protein-coding genes. 
Next, we use divergence dating to statistically test the 
temporal prediction of both the MR and the WC hy- 
potheses that recent (Pleistocene or possibly Pliocene) 
climatic cycles are a major force promoting Brookesia 
species diversification. 

Major climatic cycles have of course occurred repeat- 
edly throughout the earth’s history, and each of these 
3 general hypotheses make spatial predictions that can 
be tested independently of any temporal predictions 
(e.g., Miocene WCs could generate the same general 
fragmentation patterns as the proposed Pleistocene con- 
tractions). Specifically, the Wilmé et al. (2006) WC model 
states that contraction of mesic habitats within adjacent 
lowland watersheds during periods of glaciation cre- 
ated gaps in ancestral species’ distributions, leading to 
allopatric differentiation/speciation. Assuming stability 
in the relative geographic positions of populations, this 
model therefore predicts that sister species should gen- 
erally occupy adjacent COEs or possibly a COE and an 
adjacent RDW (although the species in the RDW would 
be expected to have a wider geographic range). Because 
most of the major watersheds are composed of several 
smaller drainages, this model could also be construed 
to predict that sister-species pairs would occupy the 
same drainage. Likewise, the MR model of Raxworthy 
and Nussbaum (1995) predicts that sister taxa will tend 
to occupy montane forested areas separated by lower 
altitude dry forests. Finally, the RB model (Pastorini 
et al. 2003) predicts species’ ranges to be delimited by 
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major lowland rivers. To evaluate these hypotheses, we 
first reconstruct species ranges using a georeferenced 
database of Brookesia distribution records and adjust 
these ranges (when sample sizes permit) through en- 
vironmental niche modeling. We then use simulations 
to test the fit of our data to the spatial predictions of 
the WC hypothesis, and we assess the qualitative fit of 
our data to the MR and RB hypotheses (all in a phylo- 
genetic context). 

Finally, each of the 3 hypotheses suggests predictions 
about the altitudinal distribution of species richness. 
The WC and RB models both predict a greater number 
of lowland species; low-headwater watersheds are 
more likely than high-headwater watersheds to become 
fragmented from their neighbors during periods of in- 
creased aridity, and rivers are widest at lower elevations 
and therefore more likely to constitute significant bar- 
riers to gene flow as altitude decreases. Conversely, the 
MR model predicts higher species diversity in mon- 
tane areas because it is here that populations tended 
to become isolated during periods of forest contrac- 
tion; lowlands were subsequently recolonized by some 
species. We use these predictions to evaluate relative 
support for these hypotheses in a nonphylogenetic con- 
text by correlation of species richness with altitudinal 
range. 


MATERIALS AND METHODS 
Taxa and Gene Regions Sampled 


A previous phylogenetic study (Townsend and Larson 
2002) identified 9 clades with uncertain interrelation- 
ships that diverged from each other early in the his- 
tory of Chamaeleonidae. Because current taxonomy 
places 2 of these major clades within the genus Brookesia, 
outgroup samples from all 7 non-Brookesia chamaeleonid 
clades were included in this study. Sampling within 
Brookesia was broad, including more than three-quarters 
of the named species as well as several undescribed 
forms, and species were sampled from multiple locali- 
ties whenever possible. To facilitate divergence-dating 
calibrations, outgroup sampling included representa- 
tives of all major squamate lineages (Table 1). 

All specimens were sequenced for 2 nuclear (RAGI, 
~1500 bp and CMOS, ~800 bp) and 3 mitochondrial 
(ND1, ~70 bp; ND2, ~1035 bp; and ND4, ~700 bp) 
protein-coding genes; we included several ND4 se- 
quences from Raxworthy et al. (2002) that expanded our 
taxonomic (3 species) or geographic (6 species) cover- 
age. Our final Brookesia sampling covered approximately 
28 species, including all but 2 of the named species. 
Museum/collection and GenBank numbers of speci- 
mens can be found in the online Supplementary Material 
(available from http://sysbio.oxfordjournals.org). 


Molecular Data and Phylogenetic Analyses 


Genomic DNA was extracted and amplified using 
standard protocols (see Townsend et al. 2004). Poly- 


TABLE 1. Nonchamaeleonid outgroup sampling for all phyloge- 
netic analyses* 


b 


Higher taxon Representative taxon 


Rhyncocephalia Sphenodon punctatus 
Squamata 
Dibamidae Dibamus 
Gekkonidae Gekko gecko 
Xantusiidae Xantusia vigilis 
Scincidae Mabuya 
Cordylidae Cordylus 
Teiidae Teiinae 
Amphisbaenia Trogonophidae 
Lacertidae Lacertidae 
Serpentes Dinodon 
Shinisauridae Shinisaurus crocodilurus 
Iguanidae Liolaemus 
Uromastycinae Uromastyx 
Leiolepidinae Leiolepis 


“All phylogenetic analyses included these outgroup taxa as well as all 
non-Brookesia chamaeleonid taxa from Figure 2. 

>’Composite taxa with more than 1 species represented among the dif- 
ferent gene partitions are called by the name of the most exclusive 
higher taxon possible. 


merase chain reaction products were sequenced. with 
ABI 3100 PRISM™ automated sequencers, and con- 
tigs were assembled using Sequencher 4.1.2 (Gene 
Codes Corporation, Ann Arbor, MI). Primer sequences 
are given in Table 2. We used the Clustal algorithm 
(Thompson et al. 1994) implemented in the program 
DAMBE (Xia and Xie 2001) to align all sequences by 
their amino acid translations and adjusted alignments 
by eye using MacClade 4.03 (Maddison D.R. and 
Maddison W.P. 2000). Ambiguously aligned regions 
were excluded from all analyses (see Results section). 
Gaps were treated as a separate (binary) character parti- 
tion in the MrBayes analyses and as missing data in the 
BEAST and maximum-likelihood (ML) analyses. 
Partitioned Bayesian and likelihood methods were 
used to infer phylogenies. BEAST (Drummond et al. 
2006) implements a Bayesian relaxed-clock method that 
allows the simultaneous estimation of topology and 
divergence times. The incorporation of a relaxed-clock 
constraint into the phylogenetic inference procedure has 
intuitive appeal, as it seems likely to model biological 
reality better than the 2 extreme alternatives of either a 
strict clock or no clock at all (i-e., the unrooted method 
of Felsenstein 1981). The relaxed clock is also expected 
to have greater statistical power due to the smaller 
number of estimated parameters relative to the un- 
rooted method (Drummond et al. 2006). MrBayes v3.1.2 
(Huelsenbeck and Ronquist 2001) and BEAST v1.4.6 
(Drummond and Rambaut 2007) were used to conduct 
Bayesian analyses under unrooted and relaxed-clock 
models, respectively. The data were divided a priori 
into partitions by gene and codon position, with mod- 
els chosen using the Akaike information criterion as 
implemented in MrModelTest (Nylander 2004). Analy- 
ses were performed under several a priori partitioning 
schemes, and the harmonic means of likelihood scores 
from the posterior distribution were compared using 
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TABLE 2. Amplification and sequencing primers used in this study 


Direction* Location Sequence Reference 

Forward ND1 5/-CAACTAATACACCTACTATGAAA-3’ Macey et al. (1997) 

Forward ND1 5’-CGATTCCGATATGACCARCT-3’ Kumazawa and Nishida (1993) 
Reverse tRNAA 5/-AAAATRTCTGRGTTGCATTCAG-3’ Macey et al. (1997) 

Forward ND4 5/-TGACTACCAAAAGCTCATGTAGAAGC-3’ Raxworthy et al. (2002) 
Reverse tRNA“ 5/-CATTACTTTTACTTGGATTTGCACCA-3’ Raxworthy et al. (2002) 
Forward RAG1 5’-TCTGAATGGAAATTCAAGCTGTT-3’ Groth and Barrowclough (1999) 
Forward RAG1 5/-CCACTTGGAAAAATACTCCCTGA-3’ This study 

Reverse RAG1 5’-GTCATCAACCAAATGTTGTATGCCTG-3’ This study 

Reverse RAG1 5'-GTGTCYACTGGGTARTCATC-3’ Townsend et al. (2004) 
Forward CMOS 5’-ATTATGCGATCMCCTMTTCC-3’ This study 

Forward CMOS 5’!-TCTGGAATTTTCTCCWTCTGT-3’ This study 

Reverse CMOS 5’-GCTACCACAGARTASAGTACA-3’ This study 


Note: tRNA = transfer RNA. 


*Mitochondrial forward and reverse primers extend the light and heavy strands, respectively. 


Bayes factors to choose a final partitioning scheme (see 
Brandley et al. 2005). 

We used the ML program RAxML (Stamatakis 2006) 
via its Web server (Stamatakis et al. 2008) to estimate a 
phylogeny and conduct nonparametric bootstrap anal- 
yses under the same partitioning scheme favored in 
the Bayesian analyses, except that all partitions were 
assigned a general time reversible model (currently the 
only option) with a proportion of invariable sites es- 
timated and branch lengths optimized separately for 
each data partition. For simplicity, further references 
to the unrooted Bayesian, relaxed-clock Bayesian, and 
maximum-likelihood analyses will be referred to by the 
initials UB, RCB, and ML, respectively. 


Divergence Time Estimates 


In most cases, fossil constraints should place the max- 
imum probability near the estimated age of the fossil, 
with probabilities dropping off rapidly for younger 
ages and more gradually for older ages (Benton and 
Donoghue 2007), effectively placing a hard bound on 
the minimum age and a soft bound (Yang and Rannala 
2006) on the maximum age. A translated lognormal 
(TL) distribution models this situation well (Hedges 
and Kumar 2004; Drummond et al. 2006; Ho 2007) and 
was used for most of the calibrations in this study. Be- 
cause overestimation of divergence times could lead to 
false rejection of Pliocene/Pleistocene speciation scenar- 
ios, we also ran analyses under a more conservative set 
of normally distributed fossil constraints in place of the 
more geologically appropriate lognormal constraints. 

Fossils or geologic data suitable for constraining 
nodes within Brookesia are not available, and we there- 
fore used several external calibration points for our 
analyses. For fossil-based calibrations, we constructed 
TL distributions in BEAST with the lower bound of the 
95% highest probability density (HPD) at a point 1% 
more recent than the estimated fossil age, and we left 
(somewhat arbitrarily) long probability tails for the soft 
maxima. The following fossil-based calibration points 
were used (Table 3): 1) the rhynchocephalian-squamate 
split within lepidosaurs (Sues and Olsen 1990; Evans 


et al. 2001), 2) the oldest known stem scincomorphs 
(here defined as the clade containing the crown taxa 
Scincidae, Cordylidae, and Xantusiidae) (Evans 1993, 
1998), 3) the oldest known anguimorph (Evans 1994, 
1998), 4) the oldest known amphisbaenian (Gao and 
Nessov 1998), and 5) the oldest teiids (Winkler et al. 
1990; Nydam and Cifelli 2002). Finally, the Comoros 
archipelago is home to the chameleon species Furcifer 
polleni and Furcifer cephalolepis, which are endemic to 
Mayotte and Grand Comoro, respectively. The oldest of 
the 2 islands (Mayotte) formed 10-15 million years ago 
(Nougier et al. 1986), and the older of these dates served 
as 6) a gamma-distributed maximum age constraint for 
the most recent common ancestor (mrca) of these taxa. 

We were concerned about the effects that substi- 
tutional saturation at deeper levels might have on 
divergence time estimates, especially with regard to 
mitochondrial DNA (mtDNA) data (Hugall et al. 2007). 
However, the mtDNA data were needed to confidently 
resolve more recent branching points (see below), sug- 
gesting also that branch lengths within Brookesia would 
be better estimated by including mtDNA data. We there- 
fore conducted analyses using 1) all data, 2) nuclear data 
plus first and second positions from the mtDNA data, 
and 3) nuclear data alone to test the robustness of our 
divergence time estimates. We also tested our fossil cal- 
ibrations using the cross-validation method of Near et 
al. (2005), which identifies potentially problematic fos- 
sils that cause incongruent age estimates of other dated 
nodes in the tree. 

All BEAST analyses were run for a sufficient number 
of generations to achieve an effective sample size of at 
least 200 for all estimated parameters, and 3 replicate 
runs were conducted for each analysis. Initial analyses 
were run without data to check the influence of the pri- 
ors on the results. BEAST output was examined for ev- 
idence of proper mixing and convergence using Tracer 
(Rambaut and Drummond 2004), runs were combined 
using LogCombiner (part of BEAST package), and max- 
imum credibility trees with divergence time means and 
95% HPDs were produced using Tree Annotator (part 
of BEAST package). All BEAST files used are included 
in the online Supplementary Material. 
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TABLE 3. Divergence date calibration priors 


Node Relevant fossil* Median (95% CI) TL zero offset 
(million years ago) (mean, SD) 
Rhyncocephalian-squamate Indeterminate taxon 231 (225-301) 224 (2.0, 1.2) 


Scincomorph-iguanian Balnealacerta 
Anguimorph-iguanian Parviraptor 
Lacertid—amphisbaenian Hodzhhakulia 
Teiid—(lacertid, amphisbaenian) Ptilotodon 
Furcifer polleni—Furcifer cephalolepis N/A 


167 (162-204) 
167 (162-204) 


161 (1.8, 1.0) 
161 (1.8, 1.0) 


102 (97-182) 97 (1.7, 1.4) 
116 (110-187) 110 (1.8, 1.3) 
6.3 (1.7-16) 0.0 (3.5, 2.0)" 


Notes: SD, standard deviation; N/A, not available. 
See text for references. 


Nonfossil calibration based on geologic data and modeled with a gamma distribution: zero offset (y-shape, y-scale) (see text). 


Distribution Mapping and Diversity Analyses 


Locality data for all Brookesia species of Madagas- 
car were gathered from museum data, our own global 
positioning system readings, and pertinent literature. 
Single-species maps that almost fully agree with the 
data used herein are reproduced in Glaw and Vences 
(2007). Small sample sizes of locality records are not 
appropriate to obtain reliable estimates of potential dis- 
tribution by modeling. For 9 species, more than 6 (range 
7-39) locality records were available. For these species, 
we defined distributions by potential distribution mod- 
els (pruned for overprediction), and for the rest of the 
species, we used point locality data. 

For the models, we used 23 variables as predictors: 
potential evapotranspiration, yearly water balance (an- 
nual rainfall minus annual evapotranspiration), number 
of months with a positive water balance (a measure of 
drought), percentage of forest cover in 2000, and 19 
climatic variables from the WorldClim database ver- 
sion 1.4 (Hijmans et al. 2005): annual mean tempera- 
ture; mean diurnal temperature range; isothermality 
(monthly /annual temperature range); temperature sea- 
sonality (standard deviation across months); maximum 
temperature of warmest month; minimum temperature 
of coldest month; annual temperature range; mean tem- 
perature of wettest, driest, warmest, and coldest quar- 
ters; annual precipitation; precipitation of wettest and 
driest months; precipitation seasonality (coefficient of 
variation); and precipitation of wettest, driest, warmest, 
and coldest quarters. 

We used Maxent v2.3 (Phillips et al. 2006) for dis- 
tribution modeling, as it performs well in predicting 
species distributions (Elith et al. 2006), even when sam- 
ple sizes are small (Hernandez et al. 2006). Analyses 
followed Vieites et al. (2008), using 1000 randomly se- 
lected data points across Madagascar as background 
pseudoabsence data. All models had areas under the 
receiver operating characteristic curve (Hanley and 
McNeil 1982) higher than 0.7, suggesting that the mod- 
els were good at discriminating between presence and 
absence sites (Fielding and Bell 1997). We applied a 
pruning algorithm to remove areas of overprediction 
from the mean model. This algorithm thresholds the 
Maxent output models using a user-defined probabil- 
ity of occurrence (t), defining a convex hull around all 
occupied regions, and buffering this hull by a given 
number of grid cells (b) and width (f), within which 


the model values are reduced until they reach zero. 
We used the following parameters: t = 40, b = 40, and 
f =80, as they provide a conservative scenario to remove 
biogeographic overprediction areas (see Kremen et al. 
2008). 

To calculate values of spatial species richness and en- 
demism, the distributional data (potential distribution 
models and point distribution data) were transformed 
into a grid cell data set and plotted on a one-quarter de- 
gree square grid covering the entire island. When more 
than 2 records were available per species, we drew a 
minimum convex polygon between localities, and grid 
cells within the polygon were considered to contain the 
species. We used Worldmap v. 4.20.18. (Williams 2002) to 
calculate species richness as the total number of species 
per grid cell. Endemism was calculated as a measure 
of range-size rarity, expressed as the percentage aggre- 
gated reciprocal range size for all species per grid cell. 


Hypothesis Testing 


Temporal concordance with the Wilmé et al. (2006) 
and the Raxworthy and Nussbaum (1995) Pliocene/ 
Pleistocene speciation scenarios (WC and MR hypothe- 
ses, respectively) was statistically tested by comparing 
the most recent tail of the 95% HPD for each sister- 
species pair in our phylogeny to the oldest limits of 
these geological epochs. The RB hypothesis was not 
formulated with an explicit time frame and thus could 
not be tested in this manner. 

The Wilmé et al. (2006) WC model predicts that sis- 
ter species should generally occupy adjacent COEs or 
possibly a COE and an adjacent RDW. To test the fit of 
our Brookesia data to the spatial aspect of this model, 
we first plotted the distributions of all Brookesia sister- 
species pairs identified in our phylogenetic analyses 
onto the Wilmé et al. (2006) watershed map (Fig. 1) 
and counted the number of pairs matching the model’s 
predictions. Next, we randomly assigned each species 
from all sister-species pairs to one of the watersheds 
collectively occupied by them and counted the number 
of pairs fitting the model’s predictions. We repeated 
this 10,000 times to generate a null distribution of ex- 
pected number of fits to the model if sister species are 
actually distributed randomly with respect to relative 
watershed positions. If less than 5% of the simulation 
rounds resulted in a number of matches equal to or 
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FIGURE 1. Map showing proposed centers of endemism (COEs; 
dark gray), retreat-dispersion watersheds (RDWs; white), and the 
3 highest peaks in Madagascar. Light gray areas are designated as 
RDWs, with possibly some function as COEs. Modified from Wilmé 
et al. (2006). 


greater than the number of matches from the real data, 
we rejected the null in favor of the Wilmé et al. (2006) 
WC hypothesis. We repeated this analysis under the in- 
terpretation that allopatric sister species occupying the 
same watershed also fit the model. 

The MR hypothesis of Raxworthy and Nussbaum 
(1995) states that rainforest contraction around high- 
elevation massifs during cooler periods isolated popu- 
lations on either side of unsuitable intervening habitat, 
leading to allopatric speciation. Referencing Figure 1, 
the 5 proposed areas of endemism (AOEs) of Raxwor- 
thy and Nussbaum (1995, see their fig. 7) correspond 
to 1) the Montagne d’Ambre mountain range in the far 
north (straddling the border between COEs 1 and 12); 2) 
a northwestern region west from Tsaratanana to Nosy 
Be and Presqu’Tle d’Ampasindava; 3) the Tsaratanana 
massif and surrounding high-altitude environs; 4) a 
northeast region extending east from Tsaratanana and 
encompassing parts of the RDWs A and a2, the southern 
tip of COE 1, and the northern section of COE 2, with 
the Antongil Bay watershed as its southwestern border; 
and finally 5) an eastern region corresponding to at least 
the southern section of COE 2. 


Several factors make it impractical to test this MR hy- 
pothesis with our data using the statistical framework 
outlined above for the WC hypothesis. First, these AOEs 
comprise only a fraction of the land area of Madagascar 
(see Fig. 1), which would necessitate the exclusion of 
several sister-species pairs that are not confined to them 
(see Results section) and thus reduce statistical power. 
Also, there are multiple potential vicariant relationships 
among most of these areas (e.g., the northeastern region 
could have once been connected by suitable habitat to 
all the other AOEs), so the number of sister-species 
pairs consistent with the MR hypothesis under the 
null model (i.e., no real effect of forest contraction on 
speciation) would be high. Furthermore, as stated in 
Raxworthy and Nussbaum (1995), there is no clear bio- 
geographic boundary between the northeastern and the 
eastern AOEs, making interpretation of relationships 
among taxa from these areas unclear. However, there 
is intuitive appeal to arguments that forest-fragment 
contraction promoted vicariant speciation in geolog- 
ically complex northern Madagascar. We thus exam- 
ine congruence of several potential examples with our 
new phylogenetic hypothesis (some of the proposed 
examples of Raxworthy and Nussbaum were flawed 
because they compared nonsister taxa; see Results 
section). 

The RB hypothesis (e.g., Pastorini et al. 2003) pre- 
dicts sister taxa to be separated by major rivers. A 
few Malagasy rivers can be unambiguously defined 
as “major” (e.g., the Betsiboka, Mangoro, and Mana- 
nara rivers) because they are quite long, with headwa- 
ters at relatively high elevations (Wilmé et al. 2006), 
and because previous works (e.g., Pastorini et al. 2003) 
have identified them as potential barriers for lemurs. 
However, for most rivers, an objective assessment of 
their potential as barriers would require combining 
parameters such as length, width, headwater eleva- 
tion, and permanence. These factors, combined with 
the difficulty of determining the size and shape of the 
test areas, preclude a formalized test of the RB hypoth- 
esis. However, we examine the spatial relationship of 
Brookesia sister taxa with several of the largest rivers to 
look for any evidence of a major role in species-level 
diversification. Our geographic sampling within some 
species also allows some evaluation of the role that 
rivers might play in intraspecific genetic structuring 
within Brookesia. 

The WC, RB, and MR hypotheses all make general 
predictions about relative species abundance at differ- 
ent altitudes. Following Wollenberg et al. (2008), we di- 
vided Madagascar into a coarser grid with cell sizes of 
82 x 63 km =5166 km? and used the ArcView extension 
“Endemicity Tools” (provided by N. Danho) to calcu- 
late Brookesia species richness and corrected weighted 
endemism (Crisp et al. 2001) for each cell. We then tested 
for nonparametric correlation of these values with al- 
titudinal range, defined as maximum-—minimum eleva- 
tions per grid cell (excluding grids with no occurrence 
of any Brookesia species). 
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RESULTS 
Data Characteristics and Phylogenetic Analyses 


Alignments were unambiguous across chameleons 
for all partitions. However, alignment of the 3’ ends 
of ND1 (~70 bp) and ND2 (~90 bp), as well as part 
of ND4 (~60 bp), was problematic across more distant 
outgroups due to high levels of amino acid replacement 
and multiple indels. We therefore coded these regions 
as missing data for nonchamaeleonid taxa. Maximum- 
likelihood analyses of near-complete ND1 and ND2 data 
across several chameleon species found very similar 
evolutionary model parameter values for these 2 genes 
(results not shown). To avoid the statistical problem of 
low data/parameter ratios, the ND1 sequence (~70 bp) 
was combined with the ND2 for partitioning purposes. 
For the same reason, data from the 3 transfer RNA genes 
separating ND1 and ND2 (and from which the loop re- 
gions were impossible to align confidently) were not 
used in this study. Bayes factor analysis on various par- 
titioning schemes (2, 3, 6, 9, and 12 partitions) strongly 
favored the 12-partition scheme (nuclear and mtDNA 
genes each partitioned by codon position; Bayes factor 
>178 in all comparisons). All results reported here are 
from this partitioning scheme. 

Within each of the ML, UB, and RCB analysis sets, 
phylogenetic results were broadly congruent across 
the different genomes and genes, with topological con- 
flicts involving only poorly supported nodes (i.e., ML 
bootstraps < 70%, Bayesian posterior probabilities [PPs] 
< 95%). One exception to this pattern involves strongly 
conflicting nuclear and mitochondrial results for the 
placement of 1 specimen of Brookesia brygooi. This dis- 
crepancy is probably due to ancient mtDNA introgres- 
sion, and the nuclear topology almost certainly reflects 
the true organismal history (see online Supplementary 
Material). Mitochondrial data for this specimen were 
therefore excluded in all further phylogenetic analy- 
ses. In general, the mtDNA data alone gave strong 
support for relationships within more terminal clades, 
but poor support for basal nodes, and the nuclear data 
showed the opposite pattern. All results presented here 
are from analyses of the combined data set (Fig. 2). 
This data set and trees from this study can be found at 
www.treebase.org under accession number SN4455. 

The mrca of the B. nasus—B. lolontany clade and all 
other Brookesia appears early in the history of the group 
(although monophyly of neither Brookesia nor this sub- 
clade is significantly supported; Fig. 2). The remainder 
of the genus is clearly monophyletic and is divided 
into 2 well-supported clades. One of these clades corre- 
sponds to the extremely miniaturized B. minima group. 
The other major clade, representing the typical Brooke- 
sia morphotype, consists of a series of 5 maximally 
supported subclades (one of which is a single species) 
whose interrelationships are strongly supported by the 
RCB analysis but less so by the UB and ML analyses 
(Fig. 2). 

Referencing the numbered clades in Figure 2, Clade 
1 contains Brookesia perarmata, Brookesia decaryi, Brookesia 


bonsi, and B. brygooi. These species are restricted to the 
dry deciduous forests of the west and/or southwest and 
together form the sister taxon of the remaining typical 
Brookesia. The eastern rainforest sister species Brookesia 
therezieni and Brookesia superciliaris (Clade 2) and the 
northern rainforest species Brookesia ebenaui (Clade 3) 
are sister taxa to progressively less inclusive clades. 
Clade 4 contains several populations that show sub- 
stantial sequence divergence from each other (11.1% 
average uncorrected ND4 divergence among specimens 
from distinct localities) but whose interrelationships are 
mostly poorly supported. The eastern rainforest species 
Brookesia thieli is paraphyletic, although basal diver- 
gences are poorly supported. One subclade of B. thieli 
appears to be the sister taxon of the morphologically 
distinctive Brookesia vadoni, suggesting the possibility of 
cryptic species within B. thieli. Finally, Clade 5 comprises 
several species (Brookesia antakarana, Brookesia ambreen- 
sis, Brookesia valerieae, Brookesia griveaudi, and Brookesia 
stumpffi) largely restricted to evergreen rainforest of 
the north, northeast, and northwest (B. stumpffi is ex- 
ceptional among Brookesia in its adaptation to both dry 
deciduous and evergreen forest). Brookesia valerieae is 
strongly supported by all analyses as the sister taxon of 
B. griveaudi. Brookesia stumpffi is weakly supported as the 
sister taxon of a strongly monophyletic B. antakarana-—B. 
ambreensis clade in the RCB analysis (Fig. 2). However, 
the other 2 methods find weak (ML, 54% bootstrap) 
to strong (UB, 96% PP) support for B. stumpffi as the 
sister taxon of the B. griveaudi—-B. valerieae clade, with 
this larger clade as the sister group to the partially 
sympatric species pair B. antakarana and B. ambreensis 
(Fig. 2). These latter 2 species exhibit very low levels 
of divergence from one another (0.6-1.3% uncorrected 
ND4 distance; see Fig. 2 inset) and may not be recip- 
rocally monophyletic. There do appear to be 2 distinct 
morphologies in this lineage (see photos in Glaw and 
Vences 2007), but mtDNA data from several additional 
specimens suggest that either incomplete lineage sort- 
ing or mitochondrial introgression is a major factor at 
this site (results not shown). 

Generally, the RCB method gave a more precise phy- 
logenetic estimate than did the UB method (83% vs. 72% 
of nodes significantly supported in Fig. 2). The increase 
in support may be a consequence of the expected in- 
crease in statistical power inherent to the relaxed-clock 
approach. Alternatively, model misspecification may 
play a role in the difference. 


Geographical Centers of Diversity and Endemism in 
Brookesia 


For the 9 modelable species, all distribution models 
predicted known localities, although comparison with 
museum records indicated some overprediction. The 
pruning algorithm removed some (but not all) areas 
appearing to have suitable environmental niches for the 
species, but which were located far from the known 
distribution ranges. This suggests that nonautecological 
(e.g., biogeographical, community ecological) factors 
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FIGURE 2. Combined mitochondrial (ND1/ND2/ND4) and nuclear (RAG1/CMOS) data, 12-partition relaxed-clock Bayesian (RCB) clado- 
gram with maximum -likelihood (ML) phylogram insert. Analyses included all ougroup taxa, although for clarity, most nonchameleon taxa are 
not shown. Major clades are color coded, and major clades of “typical” Brookesia are numbered (see text). RCB/unrooted Bayesian (UB) PP 


>95% are denoted by asterisks (PP between 90% and 95% are shown as actual values) above branches, and ML bootstrap values >50% are 
shown below branches. Black diamond indicates alternate attachment point for Brookesia stumpffi clade in the UB and ML analyses. Dashed 


terminal branches mark individuals represented by ND4 data only. 
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FIGURE 3. Maps illustrating (a) differential Brookesia species diversity (richness) across Madagascar, quantified as numbers of species present 
within one-quarter degree square grid cells covering the entire island and (b) degree of endemism in different parts of Madagascar. Endemism is 
calculated as a measure of range-size rarity, expressed as the percentage aggregated reciprocal range size for all species per one-quarter degree 


square grid cell. 


are needed to explain the restricted microendemic dis- 
tribution of some species. 

Brookesia distributions span the island, although they 
are absent from many large areas in the central high- 
lands and the south. Three main centers of diversity 
(Fig. 3a) are found in the northeast, north, and north- 
west, respectively. These areas are also COEs, with the 
highest percent values of range-size rarity (Fig. 3b). 
However, there are several areas, mainly in the west 
and northwest regions of the island, with a high de- 
gree of endemism corresponding to several species 
only known from single locality records. 

The 3 main clades recovered in our phylogenetic 
analyses (Figs. 2 and 4) show distinct biogeographic 
patterns. The basally diverging B. nasus clade contains 
2 species, B. nasus in the southeast and B. lolontany in the 
north, suggesting an old north-south connection. The B. 
minima clade shows the highest degree of endemism; 
most of the 11 species are restricted to very small areas 
in the north and north-central regions. The third and 
largest clade (typical Brookesia) spans most of the island, 
although centers of both diversity and endemism in this 
clade are in the north. Within this clade, B. brygooi, B. 
bonsi, B. decaryi, and B. perarmata (all western species; 
Clade 1, Fig. 2) form the sister taxon of the remaining 
species, none of which are restricted to the western 
forests. 


Timing of Diversification in Brookesia 


The RCB 12-partition analysis places the basal split 
within Brookesia at approximately 72 million years ago 


(95% confidence interval [CI] from ~63 to ~81 million 
years ago; Fig. 4), which is possibly older than any 
divergence within non-Brookesia chameleons (see Sup- 
plementary Material). Generally, divergences between 
recognized species groups and even sister taxa appear to 
be quite deep. Within the typical Brookesia clade, aside 
from the problematic B. ambreensis—B. antakarana com- 
plex, the most recent sister-species split (B. bonsi—B. 
decaryi) has a 95% CI extending no more recently than 
the Middle Miocene (Fig. 4). Within the B. minima clade, 
mean estimates for sister-taxon divergences are even 
older (the most recent divergence time mean is at the 
Eocene-Oligocene boundary), although the 95% CIs 
within the 2 clades overlap. Thus, we can confidently 
reject the temporal component (i.e., Pliocene—Pleistocene 
time frame) of both the MR (Raxworthy and Nussbaum 
1995) and the WC (Wilmé et al. 2006) hypotheses. 

Additional analyses using alternative calibration dis- 
tributions and data sets had little effect on divergence 
time estimates. Analyses performed with normally dis- 
tributed calibrations (with standard deviations arbitrar- 
ily set at 10% of the mean) gave average intrachameleon 
divergences about 10% more recent than those using 
lognormally distributed calibrations. Likewise, analy- 
ses that excluded mtDNA third codon positions and all 
mtDNA data gave respective divergence estimates on 
average about 15% and 18% more recent than the full- 
data analysis. In all these analyses, Brookesia sister-taxon 
divergences remained solidly Miocene in age, and thus, 
our above conclusions are unaffected. 

Using the method of Near et al. (2005), we found no 
significant incongruence among our fossil calibrations 
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FIGURE 4. Chronogram from the 12-partition (ND1/ND2, ND4, CMOS, and RAG1, each partitioned by codon position) RCB phylogenetic 
analysis. Time units on scale bar in millions of years ago. See Supplementary Materials for exact divergence time means and 95% CIs for all 
nodes. Maps above chronogram illustrate patterns of species richness (quantified as numbers of species present within one-quarter degree 
square grid cells) for the 3 main clades (1 = Brookesia nasus group, 2 = Brookesia minima group, and 3 = “typical” Brookesia group). 
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with one exception; the root calibration (representing 
the rhynchocephalian/squamate split) was a signifi- 
cant outlier (P = 0.017). Following Near et al. (2005), 
this calibration should be discarded as potentially mis- 
leading. Exclusion of this calibration point draws all 
extrachamaeleonid divergences to markedly more re- 
cent dates, and the estimated root age varies from ap- 
proximately 170-150 million years ago. However, fossil 
rhynchocephalians are known from multiple Laurasian 
and Gondwanan localities by the Late Triassic (Evans 
et al. 2001). Thus, even allowing for reasonable error in 
fossil dates, Middle to Late Jurassic estimates for the 
Rhynchocephalia-Squamata split are clearly untenable 
(see Marshall 2008 for further discussion of the poten- 
tial pitfalls of excluding outliers). We therefore favor the 
results from the full constraints analysis. That said, we 
note that even in analyses performed without the root 
constraint, in no case did the 95% HPDs for Brookesia 
sister-taxon divergences reach the Pliocene. 


Fit of Spatial Hypotheses 


We compared the distribution of Brookesia with the 
12-15 COEs recently proposed by Wilmé et al. (2006) 
(Fig. 1). Of approximately 30 Brookesia species, 18 are 
limited to single COEs (3 of these are also found in 1 
RDW), 9 are found in 2 or more separate COEs (1 is 
also found in an RDW), and 2 additional species are 
limited to single RDWs. However, these numbers are 
difficult to interpret for several reasons. First, Brookesia 
distributions are heavily skewed toward the northern 
end of the island, and certain COEs are occupied by 
disproportionate numbers of species. For example, 22 
of 30 species (73%) are distributed across just 6 COEs, 
and 7 of these species are restricted to a single COE (S. 
Bemarivo/N. Mangoro [2] on the eastern and northeast- 
ern coast) (Fig. 1). Furthermore, the level of microen- 
demism within Brookesia is extreme: 18 of 30 species 
are essentially known from single localities (Glaw and 
Vences 2007). Thus, most species of Brookesia will be 
confined to single COEs regardless of the mechanism 
of their isolation. But perhaps most importantly, sim- 
ple tabulations of species counts across proposed COEs 
ignores phylogeny completely. 

To test for spatial concordance with the WC hy- 
pothesis, we first plotted the distributions of the 10 
strongly supported (PP > 95%) Brookesia sister-species 
pairs identified in our phylogenetic analysis (i.e., 20 
species, see Fig. 2, also Supplementary Table S1) onto 
the Wilmé et al. watershed map (Fig. 1). Although 
none of these sister-species pairs are strictly confined 
to adjacent COEs, Brookesia tuberculata is found on 
Montagne d’Ambre, which straddles COEs 1 and 12, 
and its sister species Brookesia sp. “Montagne de Fran- 
cais” occupies COE 1. Two sister-species pairs (Brookesia 
exarmata [8]/ Brookesia dentata [G; also 9] and Brookesia 
sp.”Betampona” [2]/Brookesia ramanantsoai [B]) follow 
the adjacent COE/RDW pattern (Figs. 1 and 2). How- 
ever, neither of these RDW residents is a wide-ranging 
species, as the Wilmé et al. (2006) hypothesis would 


predict. Nonetheless, accepting these 3 pairs as poten- 
tial fits to the model, we performed 10,000 rounds of 
random allocations of the 20 relevant species (i.e., those 
belonging to sister-species pairs) to the watersheds 
collectively occupied by them and found that 12% of 
the simulation rounds allocated the members of 3 or 
more sister-species pairs to adjacent watersheds (i.e., 
P = 0.12). When allopatric taxa occupying the same 
watershed were included as potential matches in the 
simulations (there were none of these in the real data), 
P = 0.32. Thus, we found no support for the spatial 
component of the WC hypothesis. 

Unlike the WC hypothesis, the MR hypothesis of 
Raxworthy and Nussbaum (1995) was formulated with 
respect to 5 AOEs in the north and east of Madagascar 
(see Materials and Methods section and Fig. 1), as op- 
posed to watersheds covering the entire island. Of the 
10 Brookesia sister-species pairs from our phylogeny 
(Fig. 2), only a subset of these taxa are distributed in 
these areas and therefore useful for evaluation of this 
hypothesis. Specifically, the species pairs B. antakarana— 
B. ambreensis, B. griveaudi-B. valeriae, Brookesia karchei- 
Brookesia peyrierasi, and B. tuberculata—B. sp. “Montagne 
de Francais” are wholly confined to these AOEs. Of 
these pairs, B. griveaudi (northeast) and B. valeriae (north- 
west) both inhabit forests less than about 900 m alti- 
tude and are absent from the higher altitude (>1600 m) 
forests of the intervening Tsaratanana massif. As pro- 
posed by Raxworthy and Nussbaum (1995), these areas 
may once have been connected by an east-west corri- 
dor of even lower altitude forest that was lost during a 
period of aridification. Tsaratanana itself is home to B. 
lolontany, likely sister taxon of the only other Brookesia 
species with populations known to inhabit this alti- 
tudinal range, B. nasus of southeastern Madagascar. 
Although none of the other species pair distributions 
are consistent with this spatial model, 2 interesting 
higher level patterns within the B. minima group clade 
are potential fits. First, B. minima is restricted to the 
northwest, whereas its sister clade (B. tuberculata + 
B. sp. “Montagne des Francais”) is restricted to the 
Montagne d’Ambre region. Second, the clade formed by 
these 3 species is found only to the west of Tsaratanana, 
whereas its sister clade (B. cf. karchei + B. peyrierasi) is 
found only to the east of Tsaratanana. 

Finally, a comparison of Brookesia distribution and 
phylogeny with major rivers yields no clear support 
for the RB hypothesis. The Betsiboka River in western 
Madagascar, which has been shown to be a major barrier 
for several lemur populations (Pastorini et al. 2003), may 
possibly lay between the sister species B. exarmata (to 
the southwest) and B. dentata (to the northeast). How- 
ever, our B. dentata is the only known specimen from 
Ankaranfantsika National Park, which actually strad- 
dles the river. At any rate, B. exarmata is confined to the 
Tsingy de Bemaraha approximately 300 km (and sev- 
eral more rivers) away to the southwest. Similarly, B. 
bonsi is found at a single locality (Namoroka) about 75 
and 100 km to the west of the large Mahavavy and 
Betsiboka rivers, respectively, and its sister species 
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B. decaryi is known only from Ankarafantsika. Once 
again, the highly restricted and disjunct distributions of 
these species make any inference of a river barrier as 
the isolating mechanism highly speculative. No other 
sister-species pair (or higher level sister-clade pair) from 
our phylogeny is divided by what might reasonably be 
considered a major river. 

River barriers can of course play a role in intraspecific 
genetic structuring (e.g., Pastorini et al. 2003; Kozak et al. 
2006; Lemmon et al. 2007), and there are a few Brookesia 
species with distributions that span large rivers (B. na- 
sus, B. superciliaris, and B. thieli in the east; B. stumpffi 
and possibly B. ebenaui in the north; and B. brygooi in 
the west). To fully study the phylogeographic effects of 
these rivers will require much more thorough sam- 
pling. However, our sampling within some species does 
allow at least some preliminary comments. The western 
B. brygooi has an extensive distribution, including pop- 
ulations on both sides of the Betsiboka River. However, 
the deepest split within B. brygooi separates a population 
to the south of the river from multiple populations span- 
ning both sides of the river. Likewise, the distribution of 
B. stumpffi spans the Sambirano, Mahavavy, and Manan- 
jeba rivers (among others) in northwestern Madagascar. 
Although the deepest split does separate a popula- 
tion to the south of the Sambirano River from several 
more northern populations, relationships among the 
remaining populations do not support the Mahavavy 
and Mananjeba rivers as barriers to gene flow. Finally, 
B. nasus is found both to the north and to the south of 
the large Mananara River in southeastern Madagascar. 
However, an individual from the northernmost pop- 
ulation is nested within a clade of individuals from 
populations to the south of the potential barrier. See 
Supplementary Material for figures illustrating these 
scenarios. Thus, although further sampling could reveal 
a more complex pattern (possibly also involving smaller 
rivers), our data do not suggest a major river effect on 
genetic structuring within any of these species. 

Species richness of Brookesia was negatively correlated 
with the minimum altitude (P < 0.05) and positively 
correlated with the maximum altitude (P < 0.001) in 
each grid cell. Species richness and corrected weighted 
endemism of Brookesia furthermore showed a signif- 
icant correlation with altitudinal range per grid cell 
(P < 0.001 and P < 0.05, respectively). These results are 
not congruent with the RB and WC hypotheses, both of 
which predict a higher species richness and endemism 
at lower elevations where rivers are widest and where 
river basins forming COEs have their headwaters. In 
contrast, the results are consistent with the MR hypoth- 
esis, which predicts more species will be restricted to 
montane regions where there is a higher probability of 
isolation during dry periods. 


DISCUSSION 


Our ability to evaluate past processes that produced 
current patterns of organismal distribution and diver- 


sity is constrained and biased by the relative availabi- 
lity of climatological and geological data from different 
periods in the earth’s history. In general, Pleistocene 
climatic changes (specifically, glaciation cycles that have 
predictable effects on sea level, aridity, and tempera- 
ture) are much better understood than those from 
the Pliocene, Miocene, and earlier (Wells 2003). This 
may explain their prominence in many glaciation- 
related diversification theories. Results from phylo- 
geographic studies have challenged the importance of 
the most recent Pleistocene and Holocene glaciations 
as major drivers of species-level diversification for sev- 
eral groups (Klicka and Zink 1997; e.g., Avise et al. 
1998; Rull 2008). In Madagascar, the existence of Pleis- 
tocene glacial activity and dynamic vegetational shifts 
has been documented (Burney 1995; Vidal-Romani et al. 
2002), but well-studied examples indicating the depen- 
dence of species diversification on these processes are 
lacking. 

Most phylogeny-based biogeographical studies of 
Madagascar have centered on the relative importance of 
vicariance and dispersal as the source of Madagascar’s 
many endemic higher taxa (Yoder and Nowak 2006 
and references therein). Given the long-standing inter- 
est in Madagascar’s high level of micorendemicity, it 
is surprising that only a few phylogenetic studies have 
explicitly estimated sibling-species divergence times 
(e.g., Vences et al. 2002; Raxworthy et al. 2008; Wirta et 
al. 2008). To our knowledge, none of these studies have 
found species-level divergence times compatible with 
Pleistocene—Pliocene radiations. Other mtDNA studies 
of various groups including geckos, skinks, spiders, and 
lemurs (Yoder et al. 2000; Vences et al. 2002; Schmitz et 
al. 2005; Wood et al. 2007; Jackman et al. 2008) have re- 
ported genetic divergences that are likewise incompat- 
ible with widespread Pleistocene—Pliocene speciation, 
using standard rates of mtDNA sequence evolution (but 
see Wirta and Montreuil 2008 for a possible Pliocene 
radiation of beetles). Thus, although well-sampled phy- 
logenies of many more groups are needed, existing 
evidence (including this study) suggests that recent cli- 
matological cycles have had relatively minor effects on 
levels of species richness and associated microendemic- 
ity on Madagascar. 

Theoretical considerations (e.g., Hewitt 1996) and 
empirical examples from other groups have suggested 
the importance of earlier glaciation cycles on the North 
American and European flora and fauna (Veith et al. 
2003; Ayoub and Riechert 2004; Hewitt 2004; Kadereit 
et al. 2004; Turgeon et al. 2005). Unfortunately, histori- 
cal species distributions in Madagascar have been very 
difficult to determine due to the paucity of Tertiary fos- 
sil deposits earlier than the Quaternary (Wells 2003). 
This has in turn greatly limited the ability to infer pa- 
leoclimates and historical vegetation cover that might 
help explain existing complex microendemic distribu- 
tion patterns. Nonetheless, assuming relative stability of 
drainage systems and topographical features, resultant 
phylogenetic patterns could be very similar between 
earlier and more recent climatic cycles. An evaluation 
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of the competing hypotheses thus depends largely on 
our knowledge of current distribution, habitat, and 


phylogeny. 


Watershed Contractions, Mountain Refuges, and 
Riverine Barriers 


Our test of the spatial component of the Wilmé et al. 
WC hypothesis is potentially biased in that phylogenet- 
ically closer species might be expected to have distribu- 
tions geographically closer than more distantly related 
species. This could result in a nonrandom association 
between sister taxa and adjacent watersheds, even if WC 
was not the cause of their separation. Thus, our failure 
to reject a random watershed /sister-taxon association, 
combined with a negative correlation between altitude 
and species richness, strongly suggests that this mecha- 
nism cannot explain current Brookesia distributions. On 
the other hand, because species-level divergences in this 
group are relatively deep (Late Eocene to Late Miocene), 
there is the possibility that extinctions, climate /habitat 
changes, and range expansions have obscured patterns 
that would support this hypothesis. These are problems 
inherent to the study of any older radiation (Losos and 
Glor 2003). Unless new relevant fossils (including those 
indicative of past climates) are discovered, the first 2 
problems are largely unavoidable. However, Brookesia 
may be less prone than many groups to the last poten- 
tial problem due to the probable low vagility of these 
leaf litter inhabitants. 

The Raxworthy and Nussbaum (1995) MR model, 
although restricted to northern Madagascar, seems to 
fit to our data reasonably well. Our study reveals that 
some of the sister-taxon pairs originally cited in favor 
of this model (B. ambreensis/B. valerieae and B. stumpffi 
/B. griveaudi) are actually not sister taxa but only more 
distantly related. Nonetheless, the Tsaratanana massif 
and associated montane areas do separate 2 sister taxa 
as well as 2 multispecies clades, and another clade is 
split between the Northwest and Montagne d’Ambre 
regions. We noted earlier that for any given forest frag- 
ment, multiple potential vicariance relationships often 
exist. This fact may account for some of the apparent 
fit to this model. However, independent support for 
this model is provided by the significant positive cor- 
relation of altitude and altitudinal range with species 
richness and endemism, which is the pattern expected if 
most speciation occurs by fragmentation and isolation 
in montane forest fragments. 

Rivers appear to have played a substantial role in the 
diversification of several lemur groups (Pastorini et al. 
2003; Goodman and Ganzhorn 2004; Louis et al. 2005), 
some frogs (Vieites et al. 2006), and possibly even some 
chameleons of the genus Furcifer (Raselimanana and 
Rakotomamalala 2003). However, this does not seem to 
be true for leaf chameleons, at least at the species level 
and above. It seems unlikely that individual Brookesia 
disperse over large distances, and they certainly can- 
not swim across rivers. However, given their small size 


and leaf litter microhabitat, they could quite plausibly 
float across rivers on debris washed away during heavy 
rains. 


Higher Level Phylogenetic and Biogeographical Patterns 


Wells (2003) hypothesized that the northward drift 
of Madagascar from its pre-Paleocene position south of 
the subtropical arid zone precipitated the aridification of 
most or all the island, causing the loss or great reduction 
of much of the existing biodiversity. Later entry into the 
“trade wind” zone of southeasterly winds (most likely 
sometime in the Eocene) brought rains to the east coast, 
causing the formation of the eastern humid forests. 
According to Wells (2003), this increased moisture, 
combined with flow of the warm southern equatorial 
current through the widening Mozambique Channel, 
also led to the formation of the western deciduous 
forests. Correlating this scenario with our dated phy- 
logeny allows some speculation on drivers of diversifi- 
cation and distribution at higher levels within Brookesia. 
Specifically, aridification may have fragmented the an- 
cient B. nasus clade into disjunct high montane areas, 
with expansion of B. nasus into lower altitudes of the 
eastern rainforest during the Eocene—Oligocene. Most 
diversification within both the B. minima and the typ- 
ical Brookesia clades occurred during this same period 
and was thus likely tied to the same expansion of mesic 
habitats. The striking size disparity between members 
of these 2 clades, combined. with largely overlapping 
distributions, suggests some sort of niche partitioning 
(most likely prey size and/or microhabitat usage), but 
data supporting or refuting this hypothesis are lacking. 

Our geographical analyses confirm that northern 
Madagascar is the center of species richness and en- 
demism of Brookesia (Fig. 3), as postulated by Raxworthy 
and Nussbaum (1995). However, within the clade of typ- 
ical Brookesia, the most basal splits do not divide taxa 
restricted to the north. Instead, the first 2 splits sepa- 
rate a western and eastern clade, respectively, from the 
remainder, which itself contains mostly northern and a 
few more nested eastern species. This pattern suggests 
ancient east-west dispersal or vicariance with subse- 
quent diversification within these areas and a later di- 
versification in the north. Although the B. minima group 
is speciose, it is split at its base into a north-restricted 
clade and a more southern clade. Thus, hypotheses of 
northern or southern origins for this group are equally 
parsimonious. However, within the southern clade, 
there is a clear basal split between 1 western (dry de- 
ciduous) and 2 eastern (rainforest) groups. Interestingly, 
the repeated east-west splits within Brookesia represent 
a type of within-landmass biome shift in one of the lin- 
eages, a phenomenon that was found to be quite rare in 
a recent global-wide study of plants (Crisp et al. 2009). 

In summary, higher level diversification within Brooke- 
sia was quite possibly tied to an Eocene—Oligocene shift 
to more mesic habitats on Madagascar, and our analyses 
support a substantial role for rainforest fragmentation 
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(especially in the elevationally heterogeneous north) in 
species-level diversification. A true test of any model 
purported to generally explain the complex patterns of 
diversity in the Malagasy fauna will require at least a 
combination of distributional and current climatic data 
such as that provided by Wilmé et al. (2006) and a se- 
ries of well-supported phylogenies (with reasonably 
accurate estimates of divergence times) from a vari- 
ety of taxonomic groups. As paleoclimatological and 
paleovegetational reconstructions become available, 
they will also be very valuable because Malagasy 
terrestrial fossil deposits from the Tertiary are virtually 
nonexistent. 
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